Dipolar origin of the gas-liquid coexistence of the hard-core 1:1 electrolyte model 
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We present a systematic study of the effect of the ion pairing on the gas-liquid phase transition 
of hard-core 1:1 electrolyte models. We study a class of dipolar dimer models that depend on a 
parameter R c , the maximum separation between the ions that compose the dimer. This parameter 
can vary from a± that corresponds to the tightly tethered dipolar dimer model, to R c — > oo, 
that corresponds to the Stillinger-Lovett description of the free ion system. The coexistence curve 
and critical point parameters are obtained as a function of R c by grand canonical Monte Carlo 
techniques. Our results show that this dependence is smooth but non-monotonic and converges 
asymptotically towards the free ion case for relatively small values of R c . This fact allows us to 
describe the gas-liquid transition in the free ion model as a transition between two dimerized fluid 
phases. The role of the unpaired ions can be considered as a perturbation of this picture. 

PACS numbers: 64.70.Fx, 05.10.Ln, 05.70.Jk 
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I. INTRODUCTION 

In recent years there has been an increasing interest in 
phase transitions between fluid phases of different elec- 
trolyte concentrations in ionic solutions. Two different 
regimes have been identified experimentally jij. The 
"solvophobic" regime occurs for large solvent dielectric 
constants that effectively turn off Coulombic interactions. 
Consequently, solvophobic phase transitions are primar- 
ily driven by unfavourable interactions between solute 
and solvent. This behavior is well described by the usual 
theory of nonelectrolyte solutions with short range inter- 
actions, and clearly leads to Ising-like critical behavior. 
By contrast, in the "Coulombic" regime the solvent has 
a low dielectric constant and electrostatic interactions 
between the solute ions drive the phase separation. In 
this case, the phase diagrams are quite asymmetric, and 
apparently mean-field critical behavior has been claimed 
P,l|, |J] , although Ising-mean field crossover has also been 
seen within a narrower range of temperatures around the 
critical than in non-ionic fluids || |^ Q, || . Some of these 
experimental studies suggest that there exists a new char- 
acteristic iength in these systems that competes with the 
correlation length for density fluctuations |?| [|]. 

Electrolyte systems in the Coulomb regime are often 
modelfed as charged hard spheres embedded in a uniform 
dielectric continuum (primitive models). Most studied is 
the "restricted primitive model" (RPM), in which the 
ions are of equal size. A vapor-liquid phase transition at 
very low temperatures and densities was predicted the- 
oretically 25 years ago |9) and by early computer simu- 
lation studies pi . Improvement of computer simulation 
techniques has allowed an increasingly precise determina- 
tion of the coexistence parameters fill, Ifll O, |f4| uB, |f6| . 



There have afso been recent studies of primitive models 
with asymmetry in size |l7], |l8|, [fl| and charge p0[ . 
Very recent results suggest that the critical behavior of 
the RPM belongs to the Ising universality class Jf6[ . 

From a theoretical point of view, different approaches 
have been used in order to explain the vapor-liquid tran- 
sition of the primitive models. Integral equations such as 
the mean spherical approximation (MSA) H, , as well 
as the Debye-Hiickel theory |2l| and Poisson-Boltzmann 
approaches (2^] have succesfully been applied to it. For 
the RPM, the most succesful theories are the pairing the- 
ories, that consider the ionic fluid as a mixture of bound 
pairs and free ions in chemical equilibrium [ pT| p4| ] , in the 
spirit of the Bjerrum's ideas pq| . However, in all these 
theories the transition is driven by the free ions, even 
when in same cases, as in the Debye-Huckel-Bjerrum ap- 
proach, the associated pairs are the dominant species. 
Analytical p^j and computer simulation |27]] results, on 
the other hand, show that the structure of the vapor 
phase is dominated by neutral clusters, mostly dimers 
and tetramers. Computer simulations have also demon- 
strated that the phase envelope of the RPM resembles 
that of the charged hard dumbbell model 17, B3. This 
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result has been confirmed by theoretical studies []29|, |30| 
The question that the present work examines in detail 
is the influence of ionic association on the vapor-liquid 
transition of primitive model electrolytes. In contrast to 
an earlier study jl7j , which only considered tightly bound 
dimers, here we examine a range of models with varying 
values of R Cl the maximum separation between ions in 
a dimer. R c — * oo corresponds to the Stillinger-Lovett 
description |u| of the free ion system. 

The structure of this paper is as follows. In Section 
II, we examine the microscopic structure of the coexist- 
ing phases of the ionic fluid and compare the correlation 
functions to those of a tightly tethered dimer model. In 
Section HI , an exact chemical representation of the ionic 



fluid as a mixture of associated pairs and free ions is in- 
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troduced. The role of the associated pairs on the phase 
coexistence is studied in detail in Section IV and the pa- 
per closes with discussion and conclusions. 



II. MICROSCOPIC STRUCTURE OF THE 
IONIC AND TETHERED DIMER SYSTEMS AT 
COEXISTENCE 

In this Section we analyse the role of pairing on 
gas-liquid coexistence of primitive model ionic fluids. 
We have studied by computer simulation an 1:1 size- 
asymmetric primitive model, in which the ions are mod- 
elled as hard spheres of diameters c+ and c_ , and car- 
rying charges +q and —q, respectively, embedded on a 
dielectric continuum of dielectric constant D (D = 1 for 
the vacuum). The interaction potential between two ions 
separated by a distance is given by: 



Dm 



(1) 



where Uh s {Tij) is the hard-core potential that takes the 
value +oo if < i (cr, + <jj) and otherwise. 

The size asymmetry is characterized by the parameter 
A, defined as: 



(2) 



Monte Carlo simulations in the neutral grand- 
canonical ensemble were performed, characterized by a 
temperature T and the configurational chemical poten- 
tial for a pair of unlike ions /i = 2p rea i — 3ksT In ( A + A_ ) , 
with p rea i = (//+ + and A± are the thermal de 

Broglie wavelengths of the ionic species: 



A+ = 



2irm±kBT 



(3) 



As usual, cubic boxes of length L under periodic bound- 
ary conditions are used. The long-ranged character of 
the Coulombic interactions is handled by the use the 
Ewald summation tecnique with conducting boundary 
conditions, with 518 Fourier-space wavevectors and real- 
space damping parameter k = 5. The relative error due 
to the infinite sums truncation in the electrostatic en- 
ergy is less than 10 -5 for random configurations in small 
systems J3^|, and this choice has been also validated by 
direct simulations of the RPM ]lj| . In order to speed up 
the simulations, the basic steps of the simulations (in- 
sertion and deletions of pairs of unlike ions) are biased 
following flli"| . Moreover, a fine-discretization approxima- 
tion is used: the positions available to each ion are the 
sites of a simple cubic lattice of spacing a. This method- 
ology has succesfully been applied to the RPM 14 , to 



1:1 size-asymetric primitive models Jl7| and to z:l sizc- 
asymetric primitive models |l3], and allows a speedup 
relative to the continuum calculations of a factor of 100 



for small systems. The results are almost indistinguish- 
able from the continuum ones for a discretization param- 
eter C = cr±/a = 10 @ 0, where a± = | (o+ + cr_) 
is the unlike-ion collision diameter. This value of the dis- 
cretization parameter (£ = 10) was used in the present 
study. Histogram rcwcighting techniques 1 33 and mixed- 
field finite size scaling methods |34|] were used to ob- 
tain the vapor-liquid envelopes and the effective criti- 
cal points, respectively. For the sake of comparison, we 
have also studied tightly tethered dipolar dimer systems 
p7[ p8[ , consisting of N pairs of a positive and nega- 
tive ion restricted to remain at separations ad satisfying 
c± < < 1.02(T-i-. Simulation details and some prelim- 
inary results for a range of values of A are presented in 
Ref. @. 

The unlike-ion collision diameter a± provides the basic 
length scale appropriate for defining both the reduced 
temperature and reduced density via 

T* = k B TDa±/q 2 and p* = 2Nal/L 3 (4) 

where N+ = N- = N is the particle number of each ionic 
species Jll], [H], p2| . The reduced simulation box length 
is defined similarly via L* = L/a± 1 and the reduced en- 
ergies and chemical potential as U* — UDa±/q 2 and 
H* = lW<7±/q 2 . 

The value of the asymmetry parameter considered in 
this paper is A = 3. For that case, the gas-liquid co- 
existence of the ionic fluid shows a shift in both tem- 
perature and density respect to the tethered dimer fluid 
(see Fig. [l]), which is a general feature when comparing 



ionic and tethered dimer systems 1 17 . On the other hand 



the asymmetry is not high so as to favor large chain-like 
neutral clusters, as occurs for bigger values of A [jl7|, [L8| . 
These features qualify this case to be a typical example 
for moderately asymmetric 1:1 electrolytes, including the 
RPM. 

Fig. [l] makes the similarity between the phase diagram 
of the ionic and the tethered dimer fluids clear. It also 
suggests in an indirect way that the pairing plays a deci- 
sive role on the gas-liquid coexistence in the ionic fluid. 
In order to clarify such a role, we have studied the ion-ion 
radial distribution functions gij{r) corresponding to the 
ionic systems, and the corresponding ones for the teth- 
ered dimer systems. For this purpose, we have considered 
the gas and liquid states at coexistence for a temperature 
T 0.95T C , with T c the corresponding critical temper- 
ature. The ion-ion radial distribution functions of the 
ionic systems, for T* = 0.0425 and the coexisting den- 
sities p* = 0.0076(15) and p* = 0.162(3) are plotted in 
Figs, g and ^, respectively. In both cases, the unlike-ion 
radial distribution function becomes very large close to 
contact, indicating the association in bound pairs of un- 
like ions. Moreover, the like-ion radial distribution func- 
tions show maxima around r* = 1.5 (in the case of g , 

this peak coincides with the contact value), and there is a 
secondary maximum in g_| for r* w 2.5. These observa- 
tions allow us to conclude that there are high correlation 
between pairs of associated unlike ions. We recall that 
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the range of densities in which the gas-liquid coexistence 
occurs prevents packing effects, so the structure is com- 
pletely given by the Coulombic interactions. 

The ion-ion radial distribution functions for the teth- 
ered dimer fluid are plotted on the Fig. ^ for the gas 
branch (T* = 0.0405, p* = O.Of 19(15)) and Fig. |for the 
liquid branch (T* = 0.0405, p* = 0.206(2)). The com- 
parison between the ionic fluid and tethered dimer fluid 
microscopic structures confirms the qualitative similarity 
between both systems. A further test of this similarity 
is found in the comparison of the neutral cluster popula- 
tions in the gas branch. We use Gillan's definition of a 
cluster [ p6| . Two ions i and j are directly bound when 
the distance between them is less than R%. This condi- 
tion defines mathematically an equivalence relationship, 
and the equivalences classes in which the ions group are 
the clusters. As the interaction between like ions is re- 
pulsive, the dependence of the cluster definition on R'L + 
or R c __ should be very weak (if they are taken smaller 
than the mean distance between two like ions). On the 
other hand, the cluster definition is not very sensitive to 
the value of if that value lies between the first min- 
imum of g^ and the mean distance between aggregates. 



In this work we have used R+ + 



R c __ = R c + _ = 1.5cr±. 



As expected, the microscopic structure of the gas phase is 
dominated by the A- ion neutral clusters |l7], ^7). Their 
fractions /n for the ionic and tethered ion systems at 
the gas phase in coexistence for T — 0.95T C are quite 
similar (see Fig. ^|), although the tethered dimer systems 
have slightly higher fractions of neutral clusters than do 
ionic systems, specially for large N, in agreement with 



the results previously reported 17 



Despite the similarities found between the ionic fluid 
and the tethered dimer fluid, there are some differences 
that can rationalize the quantitative differences between 
their phase diagrams. The unlike- ion radial distribution 
functions of the tethered dimer fluid differ qualitatively 
from the ionic counterparts close to the contact value, 
since in the former there is a jump at the maximum al- 
lowed separation of two ions of a dimer, while the ionic 

<7_| takes smoothly higher values on that range of values 

of r* . Consequently, the condition on the confinement of 
the ions that compose a tethered dimer must be relaxed 
in order to refine the pairing concept in the ionic flu- 
ids. In the next Section we will introduce an appropriate 
framework for such a goal. 



III. CHEMICAL PICTURE OF IONIC FLUIDS 

We consider an ionic fluid (q + = — q_ = q) be con- 
tained in a volume V and in equilibrium with a reser- 
voir at a temperature T and chemical potentials /i+ and 
/!_ . For simplicity, the neutral grand-canonical ensemble, 
in which only neutral configurations are allowed, will be 
considered. This ensemble has been shown to be equiv- 
alent to the usual grand-canonical ensemble in the ther- 
modynamic limit |35|| . The neutral grand-canonical par- 



tition function can be written as: 

AT/ [w 



OO 

E 



Z v [Ni,Ni) (5) 



where iVj is the number of ions of each species, (3 = 
l/(fcsT), p = (p + + P-)/2, A + and A_ are the thermal 
de Broglie wavelengths corresponding to each species, 
and Zu(Ni, N) is the canonical configurational partition 
function corresponding to a fixed number of ions at the 
temperature T and enclosed in the volumen V: 

Zu(N t ,N t ) = f dr+drr 



■ dr Ni dT Ni 



(6) 

with U = UP ,l « s = U hc + J2i<j QiQj/{Drij) the total po- 
tential energy and Uhc the hard-core contribution, equal 
to +oo if two particles overlap, and otherwise (other 
tempered potentials can be used, but the main results of 
this Section will remain unchanged). If the fluid particles 
are strongly associated into bound (+, — ) pairs, as occurs 
in the low-temperature and low-density region in which 
the vapor-liquid transition occurs for the ionic fluids, the 
"physical" representation described above will be inad- 
equate, and it can replaced by a "chemical" picture, in 
which the fluid is composed by associated pairs and free 
ions. First, a rule that unequivocally identifies bound 
pairs for each configuration of the ionic fluid (up to a set 
of null measure in the configurational space) is needed. 
Once such a rule is defined, we can write the canonical 
configurational partition function as: 



Ni 



Z u (N i ,N i )= ]T Z^(N f ,N f ,N p ) 



(7) 



N„=0 



where N p is the number of bound pairs, Nf = Ni — 
N p is the number of free ions of each species, and 
Zy{Nf, Nf, N p ) is the canonical configurational partition 
function corresponding to a system of N p associated pairs 
and Nf free ions of each species: 

Zu(Nf,Nf, N p ) - ^TT^ .„ d (l)--- d (N p ) 



N p \ (Nf\) 2 J v 2(N p + N f ) 

x drfdr^ . . . dr+^dr^ e' 



-pu* 



(8) 



where (i) = {r + ,r }j are the positions of the ions that 
compose the associated pair i, and r^ - '' corresponds to 
the coordinates of a +(— ) free ion i, respectively. The po- 
tential energy U* = U chem does not coincide, in general, 
with the physical potential energy U phys , since different 
configurations of the "chemical" mixture can be compat- 
ible with a given ionic configuration. 

Substituting equation (f7|) into the equation (Q) and 
rearranging the resulting expression, the neutral grand- 
canonical partition function can be written as: 



N f =0N p =0 v + 7 



E 

N p =l 

x ZZ(N f ,Nf,N p 



A^ J 



p 



(9) 
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with fi p = 2/i and A p — ^/A + A_. These last defini- 
tions correspond to the chemical equilibrium conditions 
between the free ions and the associated pairs |36| . 
We remark that this derivation is independent of the cri- 
terion chosen to define the pairs. However, the choice 
must be such that matches the microscopic structure of 
the fluid. As we have seen in the previous Section, the va- 
por structure is dominated by neutral aggregates, mostly 
dimers and tetramers. This fact is consistent with the 
Bjerrum ideas of pairing pq] . Consequently, a distance- 
based criterion seems to be the most appropriate: two 
unlike ions that are closer than R c (R c being a suitable 
cutoff distance) are considered as an associated pair. It is 
easy to see that such a criterion does not define uniquely 
the associated pairs when the population of tetramers 
and higher order neutral clusters is not negligible, as it 
can be seen from Fig. [?], since there are different arrange- 
ments of associated pairs that are compatible with the 
same ionic configuration. Hence, a more systematic cri- 
terion is needed, in order to be able to define associated 
pairs for almost every ionic configuration, while keeping 

I 



where U phys is the physical potential energy due to the 
hard-core and eletrostatic interactions, and the other 
terms correspond to effective interactions of entropic ori- 
gin needed to reduce the configurational space of the 
"chemical" system. The term v p ({r + , r~}i) depends only 
on the positions of the ions that compose an associated 
pair and confines them to be at a distance shorter than 
R c - 

f0 |r+-rr|<£ c 

«*({r+ r-}<) = 4 (11) 
[+oo \r+-r;\>R c 

The pairwise potential energy v pp ({r + , r~}i, {r + , r~}j) 
corresponds to an steric hindrance condition between 
two associated pairs pl| , since the distances between 
two unlike ions corresponding to different pairs cannot 
be shorter than the minimum distance between the ions 
that compose each pair: 

f+oo \Tf-rf\<dij 
«»({r+ r-} i( {r+ r-} 3 -) = I 

I otherwise 

(12) 



the intuitive definition of an associated pair as encom- 
passing unlike ions that are at the closest distances. We 
use a suitable modification of the Stillinger-Lovett pair 
definition on a given ionic configuration |31| . In this pre- 
scription, all the distances between two unlike ions are 
computed, and then the first pair is defined as the two 
unlike ions at closest distance. This step is repeated, 
taking into account only ions that remain unpaired from 
previous steps for the evaluation of (+, — ) distances. We 
stop when the minimum distance between two unlike ions 
is greater than R c , considering the remaining ions as free 
ions (no free ions were considered in plj] , and conse- 
quently R c = oo in that case). This protocol provides 
an unique configuration of associated pairs and free ions 
for almost each ionic configuration, since the method is 
ambiguous only in a subset of the ionic configurations in 
which at least two (+,— ) distances are equal, and the 
measure of such a set is null in the configurational space. 

Obviously, different pairing prescriptions can be used. 
However, this protocol provides an explicit expression for 
the "chemical" potential energy Jj chem : 



(10) 



with dij = min(|r+ - r~|, |rt - rT|). 

The interaction between the associated pairs and the 
free ions is modified by the term u p± ({r + , r~}.;, ry), that 
prevents the free ion to be closer to the unlike ion of the 
associated pair than its partner: 

f+co \rf-rf\ <K+-rr| 
v p± ({r+,r-} l ,rf)=l 

I otherwise 

(13) 

Finally, two unlike free ions cannot be at a shorter 
distance than R c due to the term u H (r+ , rj): 

f+oo \r+-rj\<R c 
v + -(r+,rj) = l (14) 
[0 \ r f-Tj\>R c 

It is not hard to see that a mixture of associated pairs 
and free ions with a potential energy given by the Eq. 
( |j~0| ) is completely equivalent to the original ionic sys- 
tem. It is interesting to note that, except the term on 
v p that ties the ions that compose an associated pair, 
the other terms are pairwise, short-ranged modifications 



N p 

jjchem = l^»«+g w P({ r +,r-} i )+ « PP ({r + ,i-} l ,{r+,r-} J ) 

i=l l<i<j<N p 
N p N f N p N f N f Nf 

j — 1 i—1 J — 1 i—1 J — 1 

i 
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of the hard-core conditions, so the potential energy of 
an allowed mixture configuration is the same as in the 
corresponding ionic configuration. Moreover, every con- 
figuration of associated pairs and free ions obtained from 
an ionic configuration by the pairing procedure described 
above fulfil the constrains induced by the added poten- 
tials. Conversely, the mixture configuration is the same 
as the one obtained from the corresponding ionic config- 
uration by the Stillinger-Lovett protocol. 

This analysis provides an exact chemical representa- 
tion of the ionic system as a mixture of bound pairs and 
free ions in chemical equilibrium. The equivalence be- 
tween representations of the system is not only at the 
level of thermodynamic properties, but also in the mi- 
croscopic structure, unlike previous studies based on the 
matching of the "chemical" and the "physical" free ener- 
gies, e.g. via the virial coefficients 37, 3j|. In those stud- 
ies the effective potentials given by Eqs. ([ll]) and ( |l4| ) 
could be guessed, but the other terms (that arise from 
matching third and higher order virial coefficients) were 
not clearly identified. Furthermore, they have not being 
used at all in pairing theories of electrolytes. We must 
stress the importance of the effective new term given by 
Eq. ( p"2| ) in the potential energy, that reduces the con- 
figuration space available to the associated pairs. This 
effect is specially important for high values of R c . An 
extreme case case that illustrates the effect of missing the 
steric hindrance term is the loosely tethered dimer fluid, 
in which the ions that compose the pair can be at any rel- 
ative distance greater than a± . In this case, the canonical 
partition function of the dimer fluid is N p \ x Zjj(N p , N p ), 
where N„ is the number of dimers, Zjj(N p , N p ) is given 
by Eq. IM) and the N p l factor is the number of different 
ways of pairing the unlike ions in each ionic configuration 
and lead to different dimer configurations. As a conse- 
quence, in the thermodynamic limit the Hclmholtz free 
energy per particle fd of the loosely tethered dimer fluid 
diverges as: 



fd 



k B T 



\nN n 



(2k 

\k B T 



1 



(15) 



where fi is the ionic Helmholtz free energy per ion, that 
it is well-behaved in the case of neutral systems |35|] . 

The association degree of the ionic fluid in the co- 
existence region has been studied by computer simu- 
lation in the framework of the previous analysis. In 
the usual grand-canonical simulations, we have identi- 
fied the associated pairs by the Stillinger-Lovett rule with 
Rc = \/SL/2, i.e. all the ions are associated, and the min- 
imum image convention is used in order to calculate the 
distances between unlike ions. The probability distribu- 
tions p(r) of having a pair an internal separation distance 
r between ions is plotted in the Fig. g for the gas and 
liquid phases at T* = 0.0443. Surprisingly, both distri- 
butions are practically identical, despite the fact that the 
corresponding densities are quite different. The distribu- 
tions show a very pronounced peak for r — a±, and a 
local maximum around r = 2.5<r±, approximately where 



presents the second local maximum. For large values 

of r, the distributions p(r) decay to zero. These results 
confirm the strong association in the ionic fluid at low 
temperatures, and that the structure is weakly affected 
by variations in density, at least in the range in which the 
gas-liquid coexistence occurs. For values of r close to the 
contact value, the distribution is well described by the 
non-interacting pair fluid probability distribution, given 
by the following expression: 



ideal 



(r) 



r 2 exp 



\Dk B Tr J 



y 2 exp 



(■B&m) d y 



1 9 



Dk B Ta± 



x exp 



Dk B Ta± 



(16) 



valid for T* <C 1 and an internal cutoff R c ^> a±. How- 
ever, for r > 1.5<7±, p(r) deviates from the ideal expres- 
sion as a consequence of the interaction between asso- 
ciated pairs. The local maximum showed by the dis- 
tribution function indicates that the most bound pairs 
are solvated by less bound pairs, to form stable neutral 
tetramers (and higher order clusters, but their inclusion 
does not seem to affect the conclusions of this Section) . 
This fact is in agreement with the features observed in 
the pair correlation functions in the previous Section. We 
must stress that the structure that p(r) presents is only 
due to the Coulombic interactions. 

It is instructive to compare our results in the coexis- 
tence region to the higher temperature ones. We have 
computed n(r) at T* = 1 and the same range of densities 
(see Fig. [|). The qualitative behavior of p(r) at high 
temperatures is similar to that predicted by Stillingcr 
and Lovett fj9| . First, it shows a maximum localized in 
r ~ On the other hand, p(r) takes significant val- 

ues in a wider range than the corresponding functions at 
lower temperatures, decaying for large values of r as r~ n , 
with n w 3. The latter prediction differs from the value 
n = 4 given in Ref. ]3l| and requires further study to 
elucidate it. 

The probability distribution p Rc (r) of having a pair an 
internal separation distance between ions when the cutoff 
distance to define a pair is equal to R c can be obtained 
from p(r) as: 



P Rc {r) 



p(r) 



P(R C ) J^ P (t)dt 



(17) 



where p(r) = p(t)dt is the cumulant distribution cor- 
responding to p(r). The density of associated pairs is 
equal to pp(R c )/2, and the total density of free ions is 
p[l — p(R c )], where p is the total density of ions in the 
physical picture. It can be seen from Fig. || that for 
R c > 3cr±, more than 95% of ions are associated into 
pairs in both vapor and liquid phases. So, if the associ- 
ated pair fluid phase separate in the same (T, p) region 
as the ionic fluid, it is expected than the free ions do play 
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a mere perturbative role in the phase coexistence. Ac- 
tually, the free ionic subsystem is a non-additive binary 
charged hard-sphere mixture, where R c plays the role of 
the unlike- ion collision diameter. However, the interac- 
tions between like ions are purely repulsive, and it is not 
expected differences between the behavior of this mixture 
and the RPM as soon as a+ and c_ are less than R c . 
Furthermore, this system is in a polar environment given 
by the associated pair subsystem, and consequently the 
effective dielectric constant D e f f of the background will 
be increased. Then the free ion subsystem is expected 
to have the same behavior as the RPM at the effective 
reduced temperature T: 



T 



D e ffk B TR c 



— — x — x T 

D o± 



(18) 



where T* is the reduced temperature. For R c > 3o±, the 
critical temperature reduces at least to one third of the 
RPM reduced critical temperature. Consequently, any 
possible free ion-driven vapor-liquid transition should 
happen at much lower temperatures, and thus the free 
ion subsystem should play no role in the vapor-liquid 
transition. In order to check this hypothesis, the phase 
diagram of the associated pair fluid has to be obtained. 
This issue will be addressed in the next Section. 



IV. PHASE BEHAVIOR OF THE ASSOCIATED 
PAIRS FLUID 

The results obtained in the previous Section suggest 
that the "chemical" picture introduced above is a very 
convenient description of the ionic fluid structure. Fur- 



J 



thermore, we can analyse the role played by the associ- 
ated pairs in the phase equilibrium by eliminating the free 
ions. We have performed grand-canonical Monte Carlo 
simulations of the associated dimer system for A = 3 and 
different values of R c . No free ions are allowed (Nf = 0) 
and the potential energy of one configuration is given by 
the Eq. ([To|). The grand-canonical free energies corre- 
sponding to these systems will consequently provide an 
upper bound to the ionic fluid grand-canonical free en- 
ergy at the same temperature and pair chemical poten- 
tial. As for the ionic system, a fine-discretization ap- 
proach with a refinement parameter £ = 10, and Ewald 
summation technique with conducting boundary condi- 
tions is used to take into account the long range charac- 
ter of the Coulombic interations. The basic steps are ei- 
ther insertion or deletion of associated pairs (chosen ran- 
domly with the same probability), biased with a Boltz- 
mann distribution that depends on the separation be- 
tween the ions that compose the dimer. An associated 
pair is inserted in the following way: the negative ion 
is placed at a random place in the box, and its coun- 
terion is placed at a relative position r± (cr± < \r±\ < 
R c ) following a probability distribution proportional to 
exp (/3o9 2 0o( r ±)) ! where — </>o( r ) is the Ewald potential 
energy between two unlike ions of unit charge. On the 
other hand, a pair is deleted with a probability propor- 
tional to exp (0oq 2 4>o( r ±)) j being r± the relative position 
of the positive ion of the pair with respect to the negative 
one. The value of /3o has been adjusted to improve the 
sampling during the simulation. In this work we have set 
0o = 8. 

The acceptance probabilities of a pair insertion 
and a pair deletion, Wfj are the following: 



W l - 

i>3 




exp(/? q 2 (/>o(r ± )) 
C 3 



exp(/3 o g 2 o (r ± )) 
C 3 



J2 exp(p q 2 Mr±)) 
fc=i 

p 

^exp(/3 O (7 2 0o(r±)) 



k=l 



(19) 



(20) 



where i and j are the initial and final configurations, re- 
spectively, N p is the number of pairs at the configuration 
i, = l/(fcsT), fi is the configurational chemical poten- 
tial, and AU is the "chemical" potential energy variation 
in the movement, including the steric hindrance terms 
given by Eq. (p"2|). As it can be seen, the acceptance 
probabilities reduce to the usual ones as 0o — > 0. 

Effective critical points for different values of L* were 
estimated by using mixed-field finite-size scaling meth- 
ods p3 and assuming Ising-like criticality. Although re- 



cent results ]4(J indicate that the pressure should also en- 
ter the field mixing, our approach should be satisfactory 
to discern the dependence of the critical parameters on 
R c . Moreover, a systematic study for the RPM case has 
shown that the assumed universality class is very likely to 
be the right one |l6| . We use histogram reweighting tech- 
niques |33| to combine the histograms from different runs 
(typically three) and estimate the critical parameters and 
their standard errors. Very long runs are needed in or- 
der to overcome the critical slowing down and the low 
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acceptance ratios due to the low temperatures involved. 
Defining a step as a try of a pair insertion/deletion, we 
have performed simulations of about 2 x 10 7 — 4 x 10 7 
equilibration steps and 2 x 10 8 — 4 x 10 8 sampling steps 
for L* = 12; 5 x 10 7 — 2 x 10 8 equilibration steps ans 
6 x 10 8 — 1.2 x 10 9 sampling steps for L* = 15; and 
2 x 10 8 equilibration steps and 1.2 x 10 9 — 1.8 x 10 9 sam- 
pling steps for L* — 18, in order to get good statistics for 
the histograms. The (effective) critical parameters are 
obtained by minimizing the deviation of the appropri- 
ately scaled mixed-field M oc p — su marginal probabil- 
ity distribution (u is the potential energy density and s is 
the mixing parameter) with respect to the corresponding 
critical 3-dimensional Ising universal function |j4|, [IT| . 
As also occurs in the ionic and tethered dimer systems, 
the matching improves as L* increases (see Fig. HO). Our 
estimations of the critical parameters for different values 
of L* and R* = R c /a± are listed in Table | and plot- 
ted in Figs, [n] and [lj. The smallest value of R* is the 
same as the maximum allowed separation in the tethered 
dimer system, and as expected the results obtained for 
the associated pairs fluid are practically indistinguishable 
from the tethered dimer system. As R* increases, it is 
observed a sharp increase in the critical temperature but 
the critical density remains almost unaffected. However, 
for R* ps 3 — 4 the critical temperature reaches a max- 
imum and then decreases, converging towards the free 
ion critical temperature. On the other hand, the crit- 
ical density also decreases towards the free ion critical 
density, although the statistical uncertainties are bigger 
for this parameter and the exact path of convergence is 
less defined. We must also point out that the statisti- 
cal uncertainties are big enough not to make possible an 
estimation of correction-to-scaling effects. 

For an heuristic explanation for such a behavior, we 
have to consider the energetically favoured configura- 
tions. It is important to consider the interaction between 
two different associated pairs, that has its minimum en- 
ergy configuration in a square ring conformation for low 
values of A, with a close-energy secondary minimum con- 
formation for the linear (+,— ) chain |PjJ. If R c is too 
small, by deforming the square it is possible to find ener- 
getically relevant configurations (i.e. the potential energy 
difference per ion of such configurations with respect to 
the minimum one is of the same order as the thermal 
energy) in which one of the associated pairs has an inter- 
nal ionic separation less than R c , but the other one does 
not fulfil such a condition. Consequently, a greater value 
of R c also increases the number of energetically relevant 
configurations between two associated pairs. Then the ef- 
fective interaction between associated pairs is enhanced, 
and thus an increase of the critical temperature should be 
expected until all the relevant configurations are allowed 
(that occurs for R c ~ 3er±). On the other hand, the av- 
erage size of the associated pairs and their aggregates is 
increased, leading to a decrease on the critical density 
(in order to keep the reduced density in terms of the ef- 
fective particle size constant). These arguments are not 



longer valid for R* > 3. Our results show the inclusion 
of associated pairs with larger internal ionic separations 
are not crucial to the gas-liquid phase transition but also 
that they bar the phase transition, as it can inferred from 
the critical parameters decrease. 

As for the ionic and tethered dimer models, histogram 
reweighting techniques |33|, [ll| allow us to obtain the 
coexistence curve up to T < 0.98T C . As the finite-size 
effects are not important far from the critical point, we 
have considered L* = 15. Taking advantage of the wide 
range of densities covered by near-critical histograms, we 
combine them with liquid subcritical state histograms in 
order to extend the density range. The extra simulations 
involve shorter runs (typically 2 x 10 8 steps after equili- 
bration). The gas- liquid coexistence curves for different 
values of R* are plotted in Fig. |l^, showing the same 
tendency as the observed one in the critical parameters. 
It is interesting to note that the vapor branches coincide 
(except close to the critical point) for R c > 3cr±. We con- 
clude from this observation that is in the liquid branch 
where the ionic fluid differs mostly from the associated 
pair system. 



V. DISCUSSION AND CONCLUSIONS 

In summary, an exact "chemical" representation of the 
ionic system as a mixture of (+, — ) associated ion pairs 
and free ions has been introduced. This representation, 
closely related to the Stillinger-Lovett pairing procedure 
for the RPM, has the advantage that not only provides 
an exact matching of between the physical and chemical 
representation thermodynamics, but also from a micro- 
scopic point of view. It also avoids an entropy catastro- 
phe that occurs for the tethered dimer model studied at 
large values of the tether length. The addition to the 
physical hamiltonian of some new pairwise hard-core in- 
teractions between the "chemical" components provides 
the suitable hamiltonian for the "chemical" representa- 
tion. In the low temperature and low density regime 
in which the ionic vapor-liquid transition occurs, such a 
representation provides a faithful characterization of the 
microscopic structure of the ionic fluid, and it can be the 
basis of new theoretical approaches. 

The analysis of the phase behavior of the system only 
composed by associated pairs indicates strongly that the 
ionic fluid vapor-liquid transition is driven by them. For 
R c < 3tr±, the critical temperature increases with R c , as 
more energetically favoured configurations are allowed. 
For R c > 3cr± , the critical parameters converge smoothly 
from above towards the ionic fluid critical parameters as 
R c increases. This fact indicates not only that the free 
ions do not drive the phase transition, but also have the 
opposite effect in the transition. The value of R c « 3a± 
corresponding to the maximum on the critical tempera- 
ture, can be regarded as the optimal size of the associated 
pair. 

Some remarks on the limitations of the present work 
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are appropriate at this point. We have focussed only on 
the effect of the associated pairs in the thermodynami- 
cal properties. However, the conducting character of the 
ionic fluid is completely driven by the free ions. On the 
other hand, the remarkable success of pairing theories 
(even when the "chemical" representation they implicitly 
use is not completely correct) remains unexplained. It is 
possible that the associated pair solvation, in addition to 
the ion-ion correlations, can mimic the pair-pair interac- 
tions. Further studies are needed to completely solve the 
origin of the ionic fluid vapor-liquid phase transition. 
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TABLE I: Dependence on R* of the estimated critical pa- 
rameters. (The 1 a statistical uncertainties refer to the last 
decimal places.) 
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FIG. 1: Gas-liquid phase diagram for the ionic system (dia- 
monds) and the tethered dimer fluid (circles) with a tether 
length equal to 1.02er±. The size asymmetry parameter is set 
A = 3. The critical parameters have been obtained from the 
extrapolation to L* — > oo of the results for L* = 12, 15 and 18 
while the subcritical coexistence curves were obtained using 
L* = 15. 




FIG. 2: Ion-ion radial distribution functions corresponding 
to the ionic system at T* = 0.0425 and p* = 0.0076(15) (the 

gas phase). In the inset, the behavior of (r*) at distances 

close to the contact value. 
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FIG. 3: The same as the Fig. 
0.0425 and p* = 0.162(3). 
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FIG. 4: Ion-ion radial distribution functions corresponding 
to the tethered dimer system at T* = 0.0405 and p* = 
0.0119(15) (the gas phase). In the inset, the behavior of 
g + -(r*) at distances close to the contact value. 
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FIG. 5: The same as the Fig. 
0.0405 and p* = 0.206(2). 



in the liquid branch: T* 



0.1- 



N 

0.01: 



A 

[] 



0.001 ^ 



2 4 6 



10 12 14 16 18 20 
N 



FIG. 6: Fraction /jv of iV-ion neutral clusters. The squares 
correspond to the ionic fluid at T* — 0.0425 and p* = 
0.0076(15), and the triangles to the tethered dimer fluid at 
T* = 0.0405 and p* = 0.0119(15). 
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FIG. 7: Different ways of pairing ions using the criterion that 
two unlike ions are paired if the distance between them is 
shorter than the cutoff R c . The cation 4 can be paired with 
either the anion 1 or the anion 2. In the first case the ions 
2 and 3 remain free, and in the latter case the ions 1 and 3 
constitute another associated pair. 
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FIG. 8: Probability distribution p(r*) of having a pair an in- 
ternal separation distance r* between ions at T* = 0.0443 and 
reduced densities p* — 0.0064(7) (solid line), corresponding to 
the gas phase, and p* = 0.149(3) (dashed line), corresponding 
to the liquid phase. The dot-dashed line corresponds to the 
ideal dimer fluid approach (see text). In the inset, the cu- 

mulant distribution p(r*) = J[ p(t)dt. The meaning of the 
symbols is the same as before. 
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FIG. 9: Probability distribution p(r*) at T* = 1. At r* = 1, 
from bottom to top the lines correspond to reduced den- 
sities p* = 0.00113(1), 0.00618(1), 0.01453(2), 0.0456(3) and 
0.1566(6). 
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FIG. 10: Rescaled marginal probability distributions p c M {x) 
for the associated pairs fluid characterized by A = 3 and 
R* — 3. The solid line is the universal function correspond- 
ing to the 3-dimensional Ising universality class, and the sym- 
bols correspond to the best-matching simulation results: the 
squares correspond to V = 12, the diamonds to L* — 15 and 
the triangles to L* = 18. 



15 



0.048 



0.046 

T * 

c 

0.044 



0.042 



10 

R * 



15 



20 



FIG. 11: Dependence of the critical temperature T* on the 
reduced cutoff distance R*, for L* — 12 (squares), L* = 15 
(diamonds) and L* — 18 (triangles). 
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FIG. 12: Dependence of the critical density p* on the reduced 
cutoff distance R*, for L* = 12 (squares), L* = 15 (diamonds) 
and L* = 18 (triangles). 
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FIG. 13: Gas-liquid coexistence curves for Rl — 1.02 (circles), 
R* — 3 (squares), R* — 5 (diamonds) and the free ion system 
(triangles) . 



